Re normalization algorithm for the calculation of spectra of interacting quantum 

systems. 
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We present an algorithm for the calculation of eigenstates with definite linear momentum in quan- 
tum lattices. Our method is related to the Density Matrix Renormalization Group, and makes use of 
the distribution of multipartite entanglement to build variational wave-functions with translational 
symmetry. Its virtues are shown in the study of bilinear-biquadratic S = 1 chains. 



PACS numbers: 75.10.Pq, 03.67.-a, 75.40.Mg 

Quantum many-body systems pose problems of great 
complexity that only in a few cases can be solved analyt- 
ically. For this reason numerical algorithms have played 
a decisive role in understanding strongly-correlated mat- 
ter. A remarkable example is the Density Matrix Renor- 
malization Group (DMRG), which is a powerful tool for 
the description of ground state properties 0. On the 
other hand, the merging of Quantum Information The- 
ory and Condensed Matter has given us a new insight 
into the physics of interacting quantum systems. The 
theory of entanglement has yielded new tools to quan- 
tify quantum correlation |3L as well as a new theoreti- 
cal framework for DMRG [4j, and helped us to develop 
algorithms to deal with problems in higher dimensions 
and also to describe time evolution,systems at finite 
temperature and quantum dissipation |a,l3- 

DMRG is a variational method over the class of Matrix 
Product States which correspond to the ID realiza- 

tion of the more general Projected Entangled-Pair States 
(PEPS) introduced in 0. In a PEPS each site in a lat- 
tice is described by a set of auxiliary systems which form 
entangled pairs with their first neighbors, and the physi- 
cal state is built by local maps onto the physical Hilbcrt 
space. The distribution of bipartite entanglement gov- 
erns the characteristics of PEPS, which are ideally suited 
to describe systems with short-range correlations, some- 
thing that explains why DMRG is particularly accurate 
in describing non-critical ground states. This observa- 
tion invites us to consider the intriguing possibility of 
modifying the auxiliary state underlying PEPS to dis- 
tribute multipartite entanglement and build new varia- 
tional classes that are more suitable for a given problem. 

In this letter we present the following results: (i) 
We define the Projected Entangled-Multipartite States 
(PEMS) and study the particular case in which a GHZ 
like state is added to the auxiliary system underlying 
PEPS. The resulting variational states have a given def- 
inite linear momentum k. (ii) The new variational class 
is used to describe efficiently excitations of translational 
invariant Hamiltonians by means of a numerical DMRG- 
likc algorithm. In this way, we can calculate the lowest 
energy branch of excitations, that is, the set of minimum 



energy eigenstates for different linear momenta {I^Jf')}- 
(iii) We also present an algorithm for the calculation of 
the sequence of excited states at a given point in momen- 
tum space, {l^l? ), 1*1 ), ■ ■ ■}, which has a broad useful- 
ness, and can also be implemented together with non- 
translational invariant DMRG-like algorithms, (iv) The 
utility of the method is shown in the study of bilinear- 
biquadratic S=l spin chains, where we find indications 
of a quantum phase characterized by nematic quasi-long 
range order, which can be realized in experiments with 
cold atoms in optical lattices licj. 

We introduce our method by considering the case of a 
chain of N d s -dimensional spins. Let us assign a set of 
auxiliary subsystems x n to each site n. In ID PEPS, only 
auxiliary systems corresponding to adjacent sites, say x n 
and x n +x, are entangled in \ ^ aU x)- To describe efficiently 
multipartite entanglement, such as the one present in 
critical or spin-wave-like states, one should consider the 
most general case in which l^aux) is multipartite entan- 
gled, that is, it cannot be reduced to a product state of 
pairs of entangled sites fll| . The physical wave-function 
is created with the aid of a product of local maps P n : 



(1) 



Each P n is a local map from x n to the spin Hilbert space 
at site n. Let us see how the choice of the proper \^ a ux) 
allows us to find a variational class with translational 
invariance. Each site n is described with the aid of two 
auxiliary subsystems, a n , b n of dimension D and a third 
subsystem c n of dimension N. a n , form a maximally 
entangled state, \<j>). On the other hand, the c n are in 
the following multipartite entangled state [l2j |: 
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n=0 



where k = rik^ir/N, with n k = 0, . . . , N — 1, and T n is 
the translation operator. The auxiliary state is, thus, 
|ijr a „ x ) = |$)®^|M fe ) (see Fig. [TJ. Note that, due to 
the multipartite entanglement in \^/ a ux), subsystems c n 
are long-range correlated, such that, after the local map- 
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FIG. 1: Schematic picture of the Projected Entangled- 
Multipartite States presented in the text. 

ping, the corresponding variational wave-function is bet- 
ter suited to describe long-range correlation. 

In this case, contrary to the non-translational ansatz 
in 0, all P n 's are given by the same operator acting 
on different sites of the chain, so that their product is 
a translational invariant operator. The local maps are 
determined by a set of Nd s matrices of dimension D, 
{Af]} a =i,...,d 3 : 

1 ' 7=1,..., N 

P = E (4r]W \s) a (a\ b (p\ c < 7 |, (3) 

such that P n acts on the auxiliary states of site n and 
returns the spin state |s) with amplitude (A? Li)a,/3, pro- 
vided a n , b n , c n are in states \a), \/3), I7), respectively. 
The resulting physical states form a variational class of 
translational invariant states with momentum k: 

N - 1 ikn 

E ^T n tr{All ] A^ ] ...A^ ] }\ Sl )...\s N ).(4) 

The physical meaning of 7 becomes clear in Eq. Q. 
Due to the presence of |M&) the auxiliary state is mapped 
onto a linear combination of ID PEPS formed by circular 
permutations of Af 1 , in which 7 represents the site in the 
chain. 

Averages of observables are computed efficiently with 
the aid of the following matrices of dimension D 2 : 

^o d = E( A N®(<-<i])*) (s,|0|s) - (5) 

s,s' 

From now on (mod TV) is implicit in all functions of 
n, d. The expectation value of any operator is a linear 
combination of TV products of D 2 x D 2 matrices: 

(O.O, ...0 N ) = ±Y1 e- ikd tv{KfKt 14 ■ ■ ■ E V' d U 6 ) 

n,d 

Eq. H15|) shows that the calculation of averages is decom- 
posed into a sum over TV Fourier components. 

Let us see how to find the state of the form J3J that 
minimizes the energy of a given Hamiltonian. We con- 
sider for concreteness the case of a short range spin 



model, H = n g M o^er^ +1 , where cr M form an orthogo- 
nal set of hermitian operators. The norm and the mean 
value of the energy can be calculated with the aid of Eq. 

m- 

(*fc|* fc ) = Y,e- ikd tr{E{> d ...E?' d } (7) 

d 

(* fc |H|*fc) = 9^- ikd ti{E 1 't.E^fE^+ 1 't.E^ d } 

Any expectation value calculated with Q is a bilinear 
form of each Afo separately, (* fc |*fc) = A]^N[n]A [n] , 
(^k\H = A[ n ^H[n]A[ n ], where A[ n ] is the vector ob- 
tained by contracting s, a, and (3 in a single index. Thus, 
the energy can be minimized with respect to the set of 
matrices A?, with a given n, by solving the generalized 
eigenvalue problem W[n]A[ n ] = eoAf[n]A[ n }. This fact al- 
lows us to find the optimum PEMS in an iterative way 
that is similar to other DMRG-like methods: once we 
have found the optimum Af^ , we replace it in the wave- 
function then we actualize 7i[n + 1], J\f[n + I], and 
repeat the minimization with respect to A^ n+1 ^ . The pro- 
cess is repeated in several sweeps from n = 1 to TV, until 
the energy converges. 

According to Eq. I|15(l TL[n] and Af[n] in i|16[) can 
be expressed as Af[n] — ^ d e~ lkd J\f[n,d], H[n) = 
Ed e_ifcrf W[n,cq, where J\f[n,d], H[n,d] depend on ma- 
trices E^ ,d . In order to find an explicit expression for 
the bilinear forms, we notice first that matrices Af, ap- 
pear twice in each product in Ea. (|15l) . both in Eq and 
EQ +d ' d , for example: 

Af[n, d] (:M ^ i} = tr{( X [ a ;,]®(Af n+d pEr 1 . d ■ *? +d ~ M x 

(Al n _ d] ® X [ a p])E? +d+1 > d ..E?- 1 > d }, (8) 

where xip] is & D x D matrix with all elements 0, but a 
1 in the entry (a, (3). A similar expression can be found 
for Ti[n, d], which includes a sum over TV terms that cor- 
responds to each interaction in the Hamiltonian, and in- 
volves also products of E™' d . 

In principle Tl[n, d] and jV[n, d] could be calculated at 
each step by means of 1)141 115|l . However the structure 
of the problem resembles that of the non-translational 
invariant calculation and is suitable for a procedure of 
initialization and actualization of block operators similar 
to the one presented in 0. First, we note that there are 
TV types of products corresponding to each Fourier com- 
ponent in 1)15(1. H[n, d] requires TV products (from the TV 
interacting terms in the Hamiltonian) of TV matrices. A 
naive estimation would thus yield TV 3 multiplications at 
each step during the optimization, but we can speed up 
the algorithm by storing products of matrices Eq" 1 in the 
proper way. For example, for building the matrix Af[l, d], 
we need the product E\ . . . E± , which is obtained 
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by calculating, and storing, a sequence of products that 
we label left operators: E x ' , E^ 1 ' , .... In the 
next step, we need E\' d . . . E d+1 ' d , and we compute it 
by multiplying the stored left operator E\' d . . . E d ' d , and 
j<jd+i,d^ wn i cn j s the fi rs t of a set of products that we 
label right operators. In the following, to build Af[n,d], 
we use the stored left operators, as well as a new right 
operator Ef +1 ' d . . . E d+n ~ 1,d , that is obtained by a sin- 
gle matrix multiplication from the right operator of the 
previous step. Finally, at n — d, we have to calculate and 
store again all the left operators, and the procedure starts 
again. A similar recipe can be used for the calculation 
of Ti[n,d\, which relies on the use of recursive relations 
between operators involving matrices E%'* 0], and the 
time for one optimization step is reduced to scale like 
iV (see the Appendix for a more detailed explanation). 
Thus, an increase of the order of N in computation time 
is required by this algorithm in comparison with the non- 
translational invariant case presented in Q . On the other 
hand if we write (0J as a Matrix Product State, we need 
D' = DN dimensional bonds, something that indicates 
that, in order to get the same accuracy than in the non- 
translational schemes, our method requires a lower D. 

In the following we show how to modify our method to 
obtain the sequence of excited states with a given linear 
momentum. Let us consider that we have previously cal- 
culated the M lowest energy eigenstates, l^^'), with lin- 
ear momentum k, described by the set of matrices B[j]r i. 
In order to proceed further and find the M + 1 lowest 
energy state, we would like to run the optimization algo- 
rithm under the following constraint: 

(*if ] |* fc ) =0, Vj = l,...,M. (9) 

Within the variational class of ID PEMS defined by Q), 
each constraint in Eq. @ is given by: 

N-l 

1*0 = E e ~ ikd ^i E m E m ■ ■ ■ ( 10 ) 

d=0 

where we use the definition: 

s 

Eq. @ implies a linear constraint in the optimization 
with respect to each Af, separately. In order to include 
this constraint in our algorithm we proceed as follows. 
At each step n we use (II II) to calculate the linear form 
n] that imposes the orthogonality to \^$) in terms 

oiA k- 

(^i**> = E ((sfc»])^mf«]W = °- (i 2 ) 

B\j, n] can also be decomposed into Fourier components, 
in the same way as Tl[n], Af[n] , and each of these com- 
ponents can be computed by the same procedure for 



actualization and storage of block operators that was 
explained below. If we contract (s, a, (3) in a sin- 
gle index, then (|12JI reads B[j, nyAi n ] = 0. The lin- 
ear constraint is incorporated to the optimization pro- 
cedure by defining projectors in the subset of states 
that are orthogonal to the M lowest states, that is, 
V[n] = l-Y, id B[iM{NB l )i,5B\j,n}\ with (N B ) itj = 
B[i,n]*B\j,n], and solving the eigenvalue problem de- 
fined by the new Hamiltonian and norm matrices given by 
H[n] -> V[n]H[n]P[n], J\f[n] -> V[n]Af[n]V[n}. The same 
idea can be used in other DMRG-related algorithms. 

We have applied our method to the study of bilinear- 
biquadratic S = 1 chains, H = J2i = 
J2i cos9SiSi+i + sm6(SiSi+i) 2 , which display a rich va- 
riety of phases depending on the parameter 9. We focus 
here on the region — 37r/4 < 9 < — 7r/2, whose char- 
acteristics have been yet not fully understood. The two 
limits, — 7r/2, — 37r/4, are exactly solvable and correspond 
to a gapped dimerized 0], and a gapless ferromagnetic 
phase |14j . respectively. So far, it has remained unclear 
what happens between the dimerized and ferromagnetic 
phases, in particular, whether the dimer order survives 
down to 9 — — 3tt/A. In Ref. ^^|, it was conjectured 
that a ID quantum nematic phase should appear as an 
intermediate phase. Since then, a few numerical works 
have dealt with this problem |l6|, |l7], |l8| , but this phase 
is yet not fully characterized. 

The ability to calculate excitations in a controlled 
way, makes our algorithm ideally suited to deal with 
this problem. Fig. |21 shows the spectrum of low energy 
states at two different points in the phase diagram. At 
9 = — 7r/2 the dispersion relation corresponds qualita- 
tively to a gapped phase. Note that for the calculation 
of the spectrum it is necessary to include the constraint 
(J5J), in order to find the second lowest energy state at 
points k = 0, k = ir. At 9 = — 0.74-7T, which lies within 
the conjectured quantum nematic phase, the spectrum 
shows a qualitative change that involves the appearance 
of a soft mode. The convergence of the algorithm with 
D is shown in Fig. [3 (a). At 6 = —ir/2 we compare 
our results for the excited state ^{z with exact results 
obtained by Bethe-ansatz . The absolute error in the 
ground state, calculated by the extrapolation of E^ to 
D = 12 (ASq ' = 2 x 10~ 3 ) agrees with the error in the 

first excited state (ASq 1 ' = 3 x 1CP 3 ) obtained from the 
comparison with the exact result at 9 = —ir/2. 

By means of Ea . (| 1 5fl we can also calculate order pa- 
rameters (OP) in the ground state. Long-range dimer 
order is characterized by a non-zero value of (T> 2 )/N 2 
in the thermodynamical limit, where T> — X]i(~l) l ^M+i 
is the bond-strength oscillation. In the interval 9 < 9 C , 
9 C « — 0.77T, our finite-size results are extrapolated to 
a value T> 2 /N 2 < 3 x 10 -5 , which is set by our accu- 
racy, estimated by comparison with higher D calcula- 
tions (Fig. El(a)). On the other hand the nematic OP is 
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-7T/2. 



given by the quadrupole tensor, whose components are 
rotations of Q zz = ^ [(Si) 2 — 2/3). Long-range ne- 
matic order is described by the isotropic squared OP, 
Q 2 = J dfl((Qff) 2 ), where Qf{ is Q zz rotated to the solid 
angle Q. We find that Q 2 (N) oc jV Q '«™ in the interval 
-3tt/4 < 6 < c with 1.4 < a„ em < 2 (see Fig. GS(b)). 
Thus, long-range nematic order is absent, in qualitative 
agreement with Coleman's theorem |l9|. The fact that 
Q 2 decays algebraically with a nem > 1 is consistent with 
the existence of quasi-long range order, as denned by al- 
gebraic decay of nematic correlation functions |2jj . Note 
that a nem evolves continuously to the value a nem = 2 in 
agreement with the exact solution at 9 = — 3vt/4 [l4| . 

The spectrum in the thermodynamical limit can be ac- 
curately determined by studying the scaling of the gaps 
with system size. In particular, the gap between k = 0, 
and k — tt is extrapolated to zero within our numerical 
accuracy (AE < 10~ 3 ) in the whole region under study 
(Fig. |U (a)). We have also studied the scaled gap be- 
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tween the lowest k — states, N {e^ — u 
quantity should grow linearly with TV in a gapped phase, 
however for values < 9 C the k = the scaled gap sat- 
urates (Fig. 21 (b))- Our results in the range < 6 C , 
9 C w — 0.77T are thus consistent either with (i) a gapless 
quantum phase with nematic quasi-long range order, (ii) 
a phase with correlation lengths longer than the size of 
the chains considered here, in which case Fig. 01 would 
not correspond to the asymptotic regime, or (iii) a gap 
that is smaller than our numerical accuracy. 

In conclusion, we have presented a DMRG-like vari- 
ational algorithm that allows us to find the lowest en- 
ergy states with a definite momentum of translational 
invariant Hamiltonians. The variational class of states 
we used in the algorithms was obtained by extending the 
concept of Matrix Product States / Projected Entangled 
Pair States to include a particular multipartite entangled 
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state. An interesting extension of this work is to explore 
how other multipartite states with long-range correla- 
tions could help in simulating e.g. critical systems. 
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Appendix: Initialization and actualization of operators. 

In this appendix we provide a few details for the efficient implementation of our numerical algorithm. We show 
that by following our procedure for initialization and actualization of operators, the time required by the algorithm in 
chains with N sites scales like N 2 . 

Let us recall the calculation of averages with the variational-wave function defined by Eq. O s ; s = (s'|0|s) are 
the matrix elements of a single site observable. We define the following matrices of size D 2 : 

E% d = J2{AU® (<-<*])*) O s , s , (13) 

s.s' 

with the following definition of tensor product: 

( A ® = A «,P B <*',0'- ( 14 ) 

The expectation value of any operator is a linear combination of N products of D 2 x D 2 matrices: 

(0,0, . . . N ) = 1 £ e-^{E n d i ElX 1 ' d ■ ■ ■ K- 14 } (15) 

n,d 

The purpose of this appendix is to show how to calculate efficiently the bilinear forms that correspond to the energy 
and norm of the variational wave-function as a function of a given Ar n -t . These bilinear forms are determined by the 
relations: 



(%i^) = ^ e - iM tr{^W d •••<•"}= E i A U)i,^L'A n ^K^^ 

d s,s' 

a,a! .fi,(3' 

<*fc|ff|**> = i £ e-* k MK l ; d EZ +1 ' d ---Er 1 - d }= E i A U)l,^ al An]{AQ a ,^ (16) 

m,d,fi s ,s' 

a,a',l3,f3' 

According to Eq. (|15|) the matrices H[n] and Af[n] in Eq. i|ltj|) can be expressed as a linear combination of Fourier 
components: 

Af[n] = E e- lkd M[n, d], H[n] = E e" iW W[n, d], (17) 

d d 

where Af[n, d], T-L[n, d] depend on matrices Eq 4 only. Remember that at step n we have to solve the generalized 
eigenvalue problem defined by W[n], jV[n], to get the optimum Af-,, and then move to A? n+1 p and calculate A/"[n], 
TL[n]. We repeat this process from n = 1 to n = N, that is, we perform a sweep. Our numerical results shows that a 
number of 10 sweeps is usually enough to converge to the minimum. 
We define the set of products: 



n.m 

e d 


= ^1 


TTtm . d 






n,m 
S n'.d 


rrin.d TTin+l.d 

— (7 


..hj x , 






,n,m 


rrin.d j-in+l.d 

= % & x 


j-,m.d 
■ ■ J 






i n.m 

V 


= Em^ 


7-in+l,d j-i7n,d | 
-fV' + • 


. . + El d . 


r?m- 



E™/) . (18) 
In terms of (|18fl we find the following expression for the norm: 
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FIG. 5: Representation of effective operators for the calculation of effective Hamiltonian and Norm operators (index d is 
omitted). 

and the effective Hamiltonian: 

WMa'/J.a'./S' = 

tr{( X $] ® (Af^)*) (Af„ +d] ® x g]) e2+^ 1 .«- 1 } + 



E^tr{(x@:]®(Af„_ d] )*j er 1 -^- 1 ^E4« +rf ]< s ®^]j ^/ +1,n - 1 }+ 

E^{(x[^]®(<- d] )j C n+< '" 1 (E^+<M®^[^]) e^ 1 '- 1 }. (20) 

In (|19l I20|l x[a] is a D x 13 matrix with all elements 0, but a 1 in the entry (a,/3). Note that we have just substituted 
A[ n ] ^ Xi wherever Af^ appears in ((191 120(1 . These expressions are equivalent to the effective Norm and Hamiltonian 
operators of DMRG. In principle, H[n, d], and JV[n, d] could be calculated at each step by means of (|14J) . This would 
imply to multiply N 3 times Eo matrices (corresponding to N sites, N interacting terms in the Hamiltonian, and N 
values of d). Since the calculation of ?i[n, d], and Af[n,d] in terms of the set of matrices A? -,, m ^ n, is the most 
time demanding part of our algorithm, it is important to find a way to reduce this number of operations. Indeed, the 
number of multiplications can be reduced to scale like N by storing and actualizing the matrices (|18() in a way that 
resembles the procedure for building block operators in DMRG. Let us see how this procedure works. 

First of all, we introduce a pictorial representation of Eqs. 1)191 120(1 . In Fig. 0we represent the operators 1(18(1 
that we need to build 7i[n,d] and J\f[n,d], at each step in the minimization of the energy, see Fig. Our task is 
to calculate efficiently these sets of products by using results of the previous steps, and recursive relations between 
operators. We notice that one needs two sets of products starting at sites n and n + d, which corresponds to each of 
the left and right terms in Fig. U3 In the following we will borrow the DMRG terminology and call these two parts 
blocks A and B. The procedure that we explain below has to be carried out independently for each of the blocks A, 
B. 

Let us consider, for example, the case of e ™ +1 > n + - 1 appearing in the A block of l(20|) . In the first step, that is, 

n = 1, we calculate and store a set of 'left operators' ' d , with n' = 2, . . . , d, something that can be done in d — 1 
steps. After solving the eigenvalue problem for im, we move to the following sites, n = 2, . . . ,d. As we calculate the 

new matrices, we actualize (but do not need to store), 'right operators' at each step n, e ^+ ld ~ n ~ 1 ; by us i n g the 'right 
operators' of the previous step n — 1. At each step 'left' and 'right operators' are combined to calculate the matrices 
that we need for (gDJl: e ™+i>»+<*-i = e n+M e d+i,n+<i-i ( gee pig The game procedure and decomposition in 'left' 

and 'right' operators has to be carried out for the products + ' n ~ , which appear in block B in Eqs. 1(191 120|) (see 
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FIG. 7: Representation of the procedure for initialization/actualization of operators in the calculation of e 
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also Fig. HJJ). 

Note that operators carry here and additional index d, corresponding to each of the Fourier components in 117fl . 
The number of operations required for actualizing j\f[n, d] at a given step is independent on N. On the other hand 
we have to calculate N contributions JV[n, d] to get Af[n] by means of Eq. (|17|) . so that the number of operations per 
step scales finally like N. On the other hand, each sweep takes N steps, so that the whole algorithm scales like N 2 . 

A procedure with the same scaling with N can be applied in the calculation of the Hamiltonian at each step n. 
We explain here, as an example, how to build the operators 1 j which appears in block A of the expression 

()2Up. In the first step, n = 1, we calculate and store all the 'left blocks' h 1 ^ ' d , e n d ' d , s" ' d ', with n' = 2, . . . , d. For this 
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FIG. 8: Calculation of the operator h n+1 ' n+d - 1 by using 'left' and 'right' operators. 



calculation we need to perform a number of multiplications proportional to d, since we can use the following recursion 
relations: e d = e d , s d ' = ii^' e d , and h d — 2^ g^£j a ^ s d + h x h d . After solving the 
eigenvalue problem for Am, we move to the following sites, n = 2, .. .,1 — d. In the following steps we actualize (but 

do not need to store), 'left blocks' at each step n, e rf +1 ' ra+ _1 , ij^ 1 '™" 1 \ h d+1 ' n+d ^ 1 , by using the 'left blocks' of the 
previous step n — 1. At each step 'left' and 'right operators' are combined to calculate the matrices that we need for 

177771, ; n+l.n+d— 1 l 7t+l,<2 d+l,n+rf— 1 . n+l,d 7d-fl,n+e£— 1 , ,n+l,d d+l,n+d— 1 /*nv 1771 

CBMd = ^ e d + e d h d s M ,rf (Fig.©. 



